## Clayton, Horrillo, and Sniderman 
## The BIAT and the AMP As Measures of Racial Prejudice in Political Science: A Methodological Assessment
## Step 3: Data analysis and visualization

# Initial settings --------------------------------------------------------

rm(list = ls())
library(tidyverse)
library(cowplot)
library(xtable)

# Read data ---------------------------------------------------------------

panel <- read.csv("output/main_panel.csv")
time <- read.csv("output/main_time.csv")


# Figure 1: Correlation between BIAT and AMP -------------------------------

# Get rid of infinite values
panel$IAT_D[panel$IAT_D=="Inf"] <- NA
panel$IAT_D[panel$IAT_D=="-Inf"] <- NA

# Destring BIAT data
panel$IAT_D <- as.numeric(panel$IAT_D)

# Run OLS regression
mod <- summary(lm(diff~IAT_D, panel))

# Create labelling shortcuts for model attributes
beta <- round(mod$coefficients[2,1], 3)
rsq <- round(mod$r.squared, 3)

# Position templates
xIntercept <- max(panel$diff)
yIntercept <- max(na.omit(panel$IAT_D))

# Text to be parsed in annotation 
beta_note <- paste("beta == ", beta)
rsq_note <- paste("italic(R)^2 == ", rsq)

# IAT vs AMP scatter
iat.amp <- ggplot(panel, aes(x = IAT_D, y = diff)) +
  geom_point(alpha = 0.1) + 
  geom_smooth(method = "lm",
              formula = "y ~ x") + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_vline(xintercept = 0, linetype = "dashed") +  
  annotate("text", x = (xIntercept + 0.235), y = yIntercept, label = beta_note, color = "black", parse = TRUE, size=6) +
  annotate("text", x = (xIntercept + 0.2), y = (yIntercept-0.11), label = rsq_note, color = "black", parse = TRUE, size=6) +
  theme_light() + 
  labs(x = "BIAT-D score", 
       y = "AMP score")

ggsave("figures/Figure1.png", iat.amp, width = 8, height = 6)
rm(list=setdiff(ls(), c("panel", "time")))

# Figure 2: Favorability --------------------------------------------------

dens.iat <- ggplot() +
  geom_density(data = panel, aes(x = IAT_D), size = 2,  color = "dodgerblue") + 
  theme_light() + 
  labs(x = "BIAT-D score (panel)", y = "Density") +
  geom_vline(xintercept = 0, size = 1.5, linetype="dashed", color = "dodgerblue")

dens.iat.2 <- ggplot_build(dens.iat)

x1 <- max(which(dens.iat.2$data[[1]]$x <= 0))

dens.iat.final <- dens.iat +
  geom_area(data=data.frame(x=dens.iat.2$data[[1]]$x[0:x1],
                            y=dens.iat.2$data[[1]]$y[0:x1]),
            aes(x=x, y=y), fill="dodgerblue", alpha = 0.5)

dens.amp.p <- ggplot(panel, aes(x = diff)) +
  geom_density(size = 2, color = "firebrick3") + 
  theme_light() + 
  labs(x = "AMP score (panel)", y = "Density") +
  geom_vline(xintercept = 0, size = 1.5, linetype="dashed", color = "firebrick3")

dens.amp.p.2 <- ggplot_build(dens.amp.p)

x1 <- max(which(dens.amp.p.2$data[[1]]$x <= 0))

dens.amp.p.final <- dens.amp.p +
  geom_area(data=data.frame(x=dens.amp.p.2$data[[1]]$x[0:x1],
                            y=dens.amp.p.2$data[[1]]$y[0:x1]),
            aes(x=x, y=y), fill="firebrick3", alpha = 0.5)

dens.amp.t <- ggplot(time, aes(x = diff)) +
  geom_density(size = 2, color = "forestgreen") + 
  theme_light() + 
  labs(x = "AMP score (time-series)", y = "Density") +
  geom_vline(xintercept = 0, size = 1.5, linetype="dashed", color = "forestgreen")

dens.amp.t.2 <- ggplot_build(dens.amp.t)

x1 <- max(which(dens.amp.t.2$data[[1]]$x <= 0))

dens.amp.t.final <- dens.amp.t +
  geom_area(data=data.frame(x=dens.amp.t.2$data[[1]]$x[0:x1],
                            y=dens.amp.t.2$data[[1]]$y[0:x1]),
            aes(x=x, y=y), fill="forestgreen", alpha = 0.5)

grid <- plot_grid(dens.iat.final, dens.amp.p.final, dens.amp.t.final, nrow = 1)
ggsave("figures/Figure2.png", grid, width = 12, height = 4)

rm(list=setdiff(ls(), c("panel", "time")))

# Table 1: Misclassification --------------------------------------------------

panel <- panel %>%
  mutate(prej.lazy = ifelse(lazydiff > 0, 1, 0),
         prej.intell = ifelse(intelldiff < 0, 1, 0),
         prej.therm = ifelse(thermdiff > 0, 1, 0),
         prej.iat = ifelse(IAT_D < 0, 1, 0),
         prej.amp = ifelse(diff < 0, 1, 0)) 

time <- time %>%
  mutate(prej.lazy = ifelse(lazydiff < 0, 1, 0),
         prej.intell = ifelse(intelldiff < 0, 1, 0),
         prej.therm = ifelse(thermdiff > 0, 1, 0),
         prej.amp = ifelse(diff < 0, 1, 0))

panel.prej.lazy <- panel %>% filter(prej.lazy == 1)
panel.prej.intell <- panel %>% filter(prej.intell == 1)
panel.prej.therm <- panel %>% filter(prej.therm == 1)
time.prej.lazy <- time %>% filter(prej.lazy == 1)
time.prej.intell <- time %>% filter(prej.intell == 1)
time.prej.therm <- time %>% filter(prej.therm == 1)

## Table 1, Row 1

tab1row1 <- cbind(prop.table(table(panel.prej.lazy$prej.iat))[2],
                  prop.table(table(panel.prej.lazy$prej.amp))[2],
                  prop.table(table(time.prej.lazy$prej.amp))[2])



## Table 1, Row 2

tab1row2 <- cbind(prop.table(table(panel.prej.intell$prej.iat))[2], 
                  prop.table(table(panel.prej.intell$prej.amp))[2],
                  prop.table(table(time.prej.intell$prej.amp))[2]) 

## Table 1, Row 3

tab1row3 <- cbind(prop.table(table(panel.prej.therm$prej.iat))[2], 
                  prop.table(table(panel.prej.therm$prej.amp))[2], 
                  prop.table(table(time.prej.therm$prej.amp))[2])

# Print and save Table 1

tab1 <- xtable(rbind(tab1row1, tab1row2, tab1row3))
print(tab1, file="tables/Table1.txt")
